# Import necessary libraries
%matplotlib notebook
get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
plt.rcParams['figure.dpi'] = 300
import glob
import pandas as pd
import numpy as np
np.random.seed(7)
from ase.build import bulk
from dscribe.descriptors import SOAP
from ase.io import read, write
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn import decomposition
from sklearn.manifold import TSNE

PATH_0 = 'POSCARS/Pz_0/*'
PATH_0p25 = 'POSCARS/Pz_0p25/*'
PATH_0p5 = 'POSCARS/Pz_0p5/*'
PATH_0p75 = 'POSCARS/Pz_0p75/*'
PATH_Ordered = 'POSCARS/Ordered/POSCAR'
OMA_0 = pd.read_csv('POSCARS/Pz_0/OMMA_Pz_0_with_Elem.txt', sep = '\t')
OMA_0p25 = pd.read_csv('POSCARS/Pz_0p25/OMMA_Pz_0p25_with_Elem.txt', sep = '\t')
OMA_0p5 = pd.read_csv('POSCARS/Pz_0p5/OMMA_Pz_0p5_with_Elem.txt', sep = '\t')
OMA_0p75 = pd.read_csv('POSCARS/Pz_0p75/OMMA_Pz_0p75_with_Elem.txt', sep = '\t')
OMA_1 = pd.read_csv('POSCARS/Ordered/OMMA_Pz_1_with_Elem.txt', sep = '\t')
OMA_last_two = pd.concat([OMA_0p75, OMA_1], ignore_index = True)
To obtain further insights into the correlation between the local atomic environment of an Fe atom and its OMA, it is desirable to have a more quantitative way, in the form of well-defined descriptors, to characterize the local atomic environments. For this purpose, we employ the smooth overlap of atomic positions (SOAP) [1] approach as implemented in the DScribe software package [2], which encodes the local atomic environment, up to a specified cutoff distance, into a rotationally invariant rep-resentation, given by the “SOAP vectors”. The length of the vectors generated by SOAP depends, among other quantities, on the number of different chemical species. In our case, we generate SOAP vectors corresponding to local atomic environments including upto first nearest-neighbors (1NN, second nearest-neighbors (2NN), and third nearest-neighbors (3NN) only. In order to visualize the correlation between the local atomic environment of the Fe atoms and their OMA, we reduce the high-dimensional SOAP vectors to two dimensions using the t-distributed stochastic neighbor embedding (t-SNE)technique [3] as implemented in scikit-learn.
#For reference, following are cutoff distances for various nearest neighbours in an FCC structure.
#rcut = 2.6 #1NN
#rcut = 3.6 #2NN
#rcut = 4.4 #3NN
#rcut = 5.1 #4NN
#rcut = 5.7 #5NN
#rcut = 6.2 #6NN
#rcut = 6.7 #7NN
#rcut = 7.2 #8NN
#rcut = 7.6 #9NN
#rcut = 8.0 #10NN
def get_SOAP_descriptor_rcut(rcut):
"""
Returns pandas dataframes (X's and Y's) where each row in X is a soap vector for an Fe atom and its corresponding value of OMA in a row in Y.
param rcut : A cutoff for local region in angstroms
"""
# Setting up the SOAP descriptor
# Check out this link (https://singroup.github.io/dscribe/latest/tutorials/descriptors/soap.html) for more information on the parameters below
periodic_soap = SOAP(
species=['Fe', 'Ni'],
rcut=rcut,
nmax=8,
lmax=6,
periodic=True,
sparse=False,
rbf = 'polynomial',
sigma = 0.5)
def soap_descriptor(PATH):
"""
Returns a list of SOAP descriptors for all selected structures POSCARS
param PATH: Path for structure files
"""
# Read directories containing structure files in ascending order
# Note: Directories are names 1, 2, 3, and so on...
all_files = glob.glob(PATH + '/POSCAR')
temp = all_files[1]
indexes = [int(temp.split(sep="/")[-2]) for temp in all_files]
order = np.argsort(indexes)
list_soap = []
for file in order:
structure = read(all_files[file])
soap_descriptor = periodic_soap.create(structure)
list_soap.append(soap_descriptor)
return list_soap
def reshape_soap_descriptor(soap_FeNi):
"""
Reshape list of SOAP descriptors from (a,b,c) --> (a*b, c)
param soap_FeNi : a list of SOAP descriptors
Returns list of soap descriptors in the form of a pandas dataframe
"""
X = np.reshape(soap_FeNi, (soap_FeNi.shape[0]* soap_FeNi.shape[1],soap_FeNi.shape[2]))
df_X = pd.DataFrame(X)
df_X = df_X.round(4)
return df_X
def process_X_Y(df_X, df_Y):
"""
Process X and Y dataframe to include unique descriptors for Fe atoms only.
param df_X: dataframe containing soap vectors for atoms
param df_Y: dataframe containing OMA values of atoms
"""
X = []
Y = []
df = pd.concat([df_X, df_Y], axis = 1)
df_Fe = df[df.Elements == 'Fe'] # Drop all Ni atoms
df_Fe = df_Fe.drop(['Elements', 'OMM_X','OMM_Z'], axis = 1) # Drop Elements column
#print("\nNumber of original rows = ", len(df_Fe))
if np.any(df_Fe.duplicated()):
print("\nIdentical rows found!\nRemoving rows\n")
df_Fe.drop_duplicates(keep = 'first', inplace = True) # Drop identical rows
#print("\nNumber of rows after removing identical data = ", len(df_Fe))
Y = df_Fe[["Z_subtract_X"]]
Y = np.array(Y)
X = df_Fe.drop(["Z_subtract_X"], axis = 1)
return X, Y
# For ordered configuration (there is only 1)
list_ordered = []
ordered = read(PATH_Ordered)
descriptor_ordered = periodic_soap.create(ordered)
list_ordered.append(descriptor_ordered)
# Create lists for all soap vectors
list_0 = soap_descriptor(PATH_0)
list_0p25 = soap_descriptor(PATH_0p25)
list_0p5 = soap_descriptor(PATH_0p5)
list_0p75 = soap_descriptor(PATH_0p75)
list_last_two = list_0p75 + list_ordered
# Convert to numpy array
soap_FeNi_0 = np.array(list_0)
#print(soap_FeNi_0.shape)
soap_FeNi_0p25 = np.array(list_0p25)
#print(soap_FeNi_0p25.shape)
soap_FeNi_0p5 = np.array(list_0p5)
#print(soap_FeNi_0p5.shape)
soap_FeNi_last_two = np.array(list_last_two)
#print(soap_FeNi_last_two.shape)
# Reshape SOAP vectors
df_X_0 = reshape_soap_descriptor(soap_FeNi_0)
df_X_0p25 = reshape_soap_descriptor(soap_FeNi_0p25)
df_X_0p5 = reshape_soap_descriptor(soap_FeNi_0p5)
df_X_last_two = reshape_soap_descriptor(soap_FeNi_last_two)
# Process
X_0,Y_0 = process_X_Y(df_X_0, OMA_0)
X_0p25,Y_0p25 = process_X_Y(df_X_0p25, OMA_0p25)
X_0p5,Y_0p5 = process_X_Y(df_X_0p5, OMA_0p5)
X_last_two,Y_last_two = process_X_Y(df_X_last_two, OMA_last_two)
return X_0,Y_0, X_0p25,Y_0p25, X_0p5,Y_0p5, X_last_two,Y_last_two
# Rcut = 2.6 (First Nearest Neighbours)
X_0_1NN,Y_0_1NN,X_0p25_1NN,Y_0p25_1NN,X_0p5_1NN,Y_0p5_1NN,X_last_two_1NN,Y_last_two_1NN = get_SOAP_descriptor_rcut(rcut = 2.6)
# Rcut = 3.6 (Second Nearest Neighbours)
X_0_2NN,Y_0_2NN,X_0p25_2NN,Y_0p25_2NN,X_0p5_2NN,Y_0p5_2NN,X_last_two_2NN,Y_last_two_2NN = get_SOAP_descriptor_rcut(rcut = 3.6)
# Rcut = 4.4 (Third Nearest Neighbours)
X_0_3NN,Y_0_3NN,X_0p25_3NN,Y_0p25_3NN,X_0p5_3NN,Y_0p5_3NN,X_last_two_3NN,Y_last_two_3NN = get_SOAP_descriptor_rcut(rcut = 4.4)
Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows Identical rows found! Removing rows
def get_t_sne_reduced_X(X, per):
return TSNE(n_components=2, perplexity=per, verbose=0, learning_rate = 300).fit_transform(X)
per = 6
# First Nearest Neighbours
X_t_sne_Pz_0_1NN = get_t_sne_reduced_X(X_0_1NN, per)
X_t_sne_Pz_0p25_1NN = get_t_sne_reduced_X(X_0p25_1NN, per)
X_t_sne_Pz_0p5_1NN = get_t_sne_reduced_X(X_0p5_1NN, per)
X_t_sne_Pz_0p75_and_1_1NN = get_t_sne_reduced_X(X_last_two_1NN, per)
# Second Nearest Neighbours
X_t_sne_Pz_0_2NN = get_t_sne_reduced_X(X_0_2NN, per)
X_t_sne_Pz_0p25_2NN = get_t_sne_reduced_X(X_0p25_2NN, per)
X_t_sne_Pz_0p5_2NN = get_t_sne_reduced_X(X_0p5_2NN, per)
X_t_sne_Pz_0p75_and_1_2NN = get_t_sne_reduced_X(X_last_two_2NN, per)
# Third Nearest Neighbours
X_t_sne_Pz_0_3NN = get_t_sne_reduced_X(X_0_3NN, per)
X_t_sne_Pz_0p25_3NN = get_t_sne_reduced_X(X_0p25_3NN,per)
X_t_sne_Pz_0p5_3NN = get_t_sne_reduced_X(X_0p5_3NN, per)
X_t_sne_Pz_0p75_and_1_3NN = get_t_sne_reduced_X(X_last_two_3NN, per)
fig, axs = plt.subplots(nrows=3, ncols=3, figsize=(12, 10),subplot_kw={'xticks': [], 'yticks': []})
# set the colormap
cmap = plt.cm.get_cmap('viridis')
# Uncomment to try different colormaps. Note: Make sure to use perceptually uniform sequential colormaps
#cmap = plt.cm.get_cmap('plasma_r')
#cmap = plt.cm.get_cmap('inferno')
#cmap = plt.cm.get_cmap('magma')
# set the marker
marker = '^'
# First nearest-neighbours
l1 = axs[0][2].scatter(X_t_sne_Pz_0p75_and_1_1NN[:,0:1], X_t_sne_Pz_0p75_and_1_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_1NN)#, cmap = cmap,edgecolors='none')
l2 = axs[0][1].scatter(X_t_sne_Pz_0p5_1NN[:,0:1], X_t_sne_Pz_0p5_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_1NN)#, cmap=cmap,edgecolors='none')
l3 = axs[0][0].scatter(X_t_sne_Pz_0_1NN[:,0:1], X_t_sne_Pz_0_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_1NN)#, cmap=cmap,edgecolors='none')
l1.set_facecolor('none')
l2.set_facecolor('none')
l3.set_facecolor('none')
# Second nearest-neighbours
l4 = axs[1][2].scatter(X_t_sne_Pz_0p75_and_1_2NN[:,0:1], X_t_sne_Pz_0p75_and_1_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_2NN)#, cmap=cmap,edgecolors='none')
l5 = axs[1][1].scatter(X_t_sne_Pz_0p5_2NN[:,0:1], X_t_sne_Pz_0p5_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_2NN)#, cmap=cmap,edgecolors='none')
l6 = axs[1][0].scatter(X_t_sne_Pz_0_2NN[:,0:1], X_t_sne_Pz_0_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_2NN)#, cmap=cmap,edgecolors='none')
l4.set_facecolor('none')
l5.set_facecolor('none')
l6.set_facecolor('none')
# Third nearest-neighbours
l7 = axs[2][2].scatter(X_t_sne_Pz_0p75_and_1_3NN[:,0:1], X_t_sne_Pz_0p75_and_1_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_3NN)#, cmap=cmap,edgecolors='none')
l8 = axs[2][1].scatter(X_t_sne_Pz_0p5_3NN[:,0:1], X_t_sne_Pz_0p5_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_3NN)#, cmap=cmap,edgecolors='none')
l9 = axs[2][0].scatter(X_t_sne_Pz_0_3NN[:,0:1], X_t_sne_Pz_0_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_3NN)#, cmap=cmap,edgecolors='none')
l7.set_facecolor('none')
l8.set_facecolor('none')
l9.set_facecolor('none')
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.88,
wspace=0.25,
hspace=0.2)
fig.text(0.65,0.9, r'$P_z \geq 0.75$',fontsize=12)
fig.text(0.42,0.9, r'$P_z = 0.5 $',fontsize=12)
fig.text(0.18,0.9, r'$P_z = 0$',fontsize=12)
fig.text(0.05,0.76, r'1NN',fontsize=12)
fig.text(0.05,0.50, r'2NN',fontsize=12)
fig.text(0.05,0.22, r'3NN',fontsize=12)
fig.text(0.07,0.86, r'(a)',fontsize=12)
fig.text(0.32,0.86, r'(b)',fontsize=12)
fig.text(0.575,0.86, r'(c)',fontsize=12)
fig.text(0.07,0.59, r'(d)',fontsize=12)
fig.text(0.32,0.59, r'(e)',fontsize=12)
fig.text(0.575,0.59, r'(f)',fontsize=12)
fig.text(0.07,0.31, r'(g)',fontsize=12)
fig.text(0.32,0.31, r'(h)',fontsize=12)
fig.text(0.575,0.31, r'(i)',fontsize=12)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
bar = fig.colorbar(l1, cax=cbar_ax, ticks=[-0.005, 0.000, 0.005, 0.010, 0.015])
bar.ax.set_yticklabels(['-0.005', '0.00', '0.005', '0.010', '0.015'], fontsize = 12) # vertically oriented colorbar
fig.text(0.99,0.54, r'OMA ($\mu_{B}$)',fontsize=12)
#plt.savefig("soap_10_unfilled.png", dpi=300, bbox_inches='tight')
plt.show()
# Let's first get the plot for P_z >= 0.75 (Fig. c)
fig, ax = plt.subplots(figsize=(4, 4))
g = ax.scatter(X_t_sne_Pz_0p75_and_1_1NN[:,0:1], X_t_sne_Pz_0p75_and_1_1NN[:,1:2], s=100, alpha=0.9, marker = '^', c = Y_last_two_1NN)
g.set_facecolor('none')
fig.colorbar(g)
# fig.savefig("soap_unfilled_0p75.png", dpi=300, bbox_inches='tight')
<matplotlib.colorbar.Colorbar at 0x7f923d9cc908>
# Set perplexity = 10
per = 10
# First Nearest Neighbours
X_t_sne_Pz_0_1NN = get_t_sne_reduced_X(X_0_1NN, per)
X_t_sne_Pz_0p25_1NN = get_t_sne_reduced_X(X_0p25_1NN, per)
X_t_sne_Pz_0p5_1NN = get_t_sne_reduced_X(X_0p5_1NN, per)
X_t_sne_Pz_0p75_and_1_1NN = get_t_sne_reduced_X(X_last_two_1NN, per)
# (Second Nearest Neighbours
X_t_sne_Pz_0_2NN = get_t_sne_reduced_X(X_0_2NN, per)
X_t_sne_Pz_0p25_2NN = get_t_sne_reduced_X(X_0p25_2NN, per)
X_t_sne_Pz_0p5_2NN = get_t_sne_reduced_X(X_0p5_2NN, per)
X_t_sne_Pz_0p75_and_1_2NN = get_t_sne_reduced_X(X_last_two_2NN, per)
# Third Nearest Neighbours)
X_t_sne_Pz_0_3NN = get_t_sne_reduced_X(X_0_3NN, per)
X_t_sne_Pz_0p25_3NN = get_t_sne_reduced_X(X_0p25_3NN, per)
X_t_sne_Pz_0p5_3NN = get_t_sne_reduced_X(X_0p5_3NN,per)
X_t_sne_Pz_0p75_and_1_3NN = get_t_sne_reduced_X(X_last_two_3NN,per)
## PLOT
fig, axs = plt.subplots(nrows=3, ncols=3, figsize=(12, 10),subplot_kw={'xticks': [], 'yticks': []})
# set the colormap
cmap = plt.cm.get_cmap('viridis')
# Uncomment to try different colormaps. Note: Make sure to use perceptually uniform sequential colormaps
#cmap = plt.cm.get_cmap('plasma_r')
#cmap = plt.cm.get_cmap('inferno')
#cmap = plt.cm.get_cmap('magma')
# set the marker
marker = '^'
# First nearest-neighbours
l1 = axs[0][2].scatter(X_t_sne_Pz_0p75_and_1_1NN[:,0:1], X_t_sne_Pz_0p75_and_1_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_1NN)#, cmap = cmap,edgecolors='none')
l2 = axs[0][1].scatter(X_t_sne_Pz_0p5_1NN[:,0:1], X_t_sne_Pz_0p5_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_1NN)#, cmap=cmap,edgecolors='none')
l3 = axs[0][0].scatter(X_t_sne_Pz_0_1NN[:,0:1], X_t_sne_Pz_0_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_1NN)#, cmap=cmap,edgecolors='none')
l1.set_facecolor('none')
l2.set_facecolor('none')
l3.set_facecolor('none')
# Second nearest-neighbours
l4 = axs[1][2].scatter(X_t_sne_Pz_0p75_and_1_2NN[:,0:1], X_t_sne_Pz_0p75_and_1_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_2NN)#, cmap=cmap,edgecolors='none')
l5 = axs[1][1].scatter(X_t_sne_Pz_0p5_2NN[:,0:1], X_t_sne_Pz_0p5_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_2NN)#, cmap=cmap,edgecolors='none')
l6 = axs[1][0].scatter(X_t_sne_Pz_0_2NN[:,0:1], X_t_sne_Pz_0_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_2NN)#, cmap=cmap,edgecolors='none')
l4.set_facecolor('none')
l5.set_facecolor('none')
l6.set_facecolor('none')
# Third nearest-neighbours
l7 = axs[2][2].scatter(X_t_sne_Pz_0p75_and_1_3NN[:,0:1], X_t_sne_Pz_0p75_and_1_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_3NN)#, cmap=cmap,edgecolors='none')
l8 = axs[2][1].scatter(X_t_sne_Pz_0p5_3NN[:,0:1], X_t_sne_Pz_0p5_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_3NN)#, cmap=cmap,edgecolors='none')
l9 = axs[2][0].scatter(X_t_sne_Pz_0_3NN[:,0:1], X_t_sne_Pz_0_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_3NN)#, cmap=cmap,edgecolors='none')
l7.set_facecolor('none')
l8.set_facecolor('none')
l9.set_facecolor('none')
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.88,
wspace=0.25,
hspace=0.2)
fig.text(0.65,0.9, r'$P_z \geq 0.75$',fontsize=12)
fig.text(0.42,0.9, r'$P_z = 0.5 $',fontsize=12)
fig.text(0.18,0.9, r'$P_z = 0$',fontsize=12)
fig.text(0.05,0.76, r'1NN',fontsize=12)
fig.text(0.05,0.50, r'2NN',fontsize=12)
fig.text(0.05,0.22, r'3NN',fontsize=12)
fig.text(0.07,0.86, r'(a)',fontsize=12)
fig.text(0.32,0.86, r'(b)',fontsize=12)
fig.text(0.575,0.86, r'(c)',fontsize=12)
fig.text(0.07,0.59, r'(d)',fontsize=12)
fig.text(0.32,0.59, r'(e)',fontsize=12)
fig.text(0.575,0.59, r'(f)',fontsize=12)
fig.text(0.07,0.31, r'(g)',fontsize=12)
fig.text(0.32,0.31, r'(h)',fontsize=12)
fig.text(0.575,0.31, r'(i)',fontsize=12)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
bar = fig.colorbar(l1, cax=cbar_ax, ticks=[-0.005, 0.000, 0.005, 0.010, 0.015])
bar.ax.set_yticklabels(['-0.005', '0.00', '0.005', '0.010', '0.015'], fontsize = 12) # vertically oriented colorbar
fig.text(0.99,0.54, r'OMA ($\mu_{B}$)',fontsize=12)
fig.suptitle('Effect of Perplexity (per = 10)', fontsize=16)
#plt.savefig("soap_10_unfilled.png", dpi=300, bbox_inches='tight')
plt.show()
# Set perplexity = 20
per = 20
# First Nearest Neighbours
X_t_sne_Pz_0_1NN = get_t_sne_reduced_X(X_0_1NN, per)
X_t_sne_Pz_0p25_1NN = get_t_sne_reduced_X(X_0p25_1NN, per)
X_t_sne_Pz_0p5_1NN = get_t_sne_reduced_X(X_0p5_1NN, per)
X_t_sne_Pz_0p75_and_1_1NN = get_t_sne_reduced_X(X_last_two_1NN, per)
# (Second Nearest Neighbours
X_t_sne_Pz_0_2NN = get_t_sne_reduced_X(X_0_2NN, per)
X_t_sne_Pz_0p25_2NN = get_t_sne_reduced_X(X_0p25_2NN, per)
X_t_sne_Pz_0p5_2NN = get_t_sne_reduced_X(X_0p5_2NN, per)
X_t_sne_Pz_0p75_and_1_2NN = get_t_sne_reduced_X(X_last_two_2NN, per)
# Third Nearest Neighbours)
X_t_sne_Pz_0_3NN = get_t_sne_reduced_X(X_0_3NN, per)
X_t_sne_Pz_0p25_3NN = get_t_sne_reduced_X(X_0p25_3NN, per)
X_t_sne_Pz_0p5_3NN = get_t_sne_reduced_X(X_0p5_3NN,per)
X_t_sne_Pz_0p75_and_1_3NN = get_t_sne_reduced_X(X_last_two_3NN,per)
## PLOT
fig, axs = plt.subplots(nrows=3, ncols=3, figsize=(12, 10),subplot_kw={'xticks': [], 'yticks': []})
# set the colormap
cmap = plt.cm.get_cmap('viridis')
# Uncomment to try different colormaps. Note: Make sure to use perceptually uniform sequential colormaps
#cmap = plt.cm.get_cmap('plasma_r')
#cmap = plt.cm.get_cmap('inferno')
#cmap = plt.cm.get_cmap('magma')
# set the marker
marker = '^'
# First nearest-neighbours
l1 = axs[0][2].scatter(X_t_sne_Pz_0p75_and_1_1NN[:,0:1], X_t_sne_Pz_0p75_and_1_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_1NN)#, cmap = cmap,edgecolors='none')
l2 = axs[0][1].scatter(X_t_sne_Pz_0p5_1NN[:,0:1], X_t_sne_Pz_0p5_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_1NN)#, cmap=cmap,edgecolors='none')
l3 = axs[0][0].scatter(X_t_sne_Pz_0_1NN[:,0:1], X_t_sne_Pz_0_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_1NN)#, cmap=cmap,edgecolors='none')
l1.set_facecolor('none')
l2.set_facecolor('none')
l3.set_facecolor('none')
# Second nearest-neighbours
l4 = axs[1][2].scatter(X_t_sne_Pz_0p75_and_1_2NN[:,0:1], X_t_sne_Pz_0p75_and_1_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_2NN)#, cmap=cmap,edgecolors='none')
l5 = axs[1][1].scatter(X_t_sne_Pz_0p5_2NN[:,0:1], X_t_sne_Pz_0p5_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_2NN)#, cmap=cmap,edgecolors='none')
l6 = axs[1][0].scatter(X_t_sne_Pz_0_2NN[:,0:1], X_t_sne_Pz_0_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_2NN)#, cmap=cmap,edgecolors='none')
l4.set_facecolor('none')
l5.set_facecolor('none')
l6.set_facecolor('none')
# Third nearest-neighbours
l7 = axs[2][2].scatter(X_t_sne_Pz_0p75_and_1_3NN[:,0:1], X_t_sne_Pz_0p75_and_1_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_3NN)#, cmap=cmap,edgecolors='none')
l8 = axs[2][1].scatter(X_t_sne_Pz_0p5_3NN[:,0:1], X_t_sne_Pz_0p5_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_3NN)#, cmap=cmap,edgecolors='none')
l9 = axs[2][0].scatter(X_t_sne_Pz_0_3NN[:,0:1], X_t_sne_Pz_0_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_3NN)#, cmap=cmap,edgecolors='none')
l7.set_facecolor('none')
l8.set_facecolor('none')
l9.set_facecolor('none')
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.88,
wspace=0.25,
hspace=0.2)
fig.text(0.65,0.9, r'$P_z \geq 0.75$',fontsize=12)
fig.text(0.42,0.9, r'$P_z = 0.5 $',fontsize=12)
fig.text(0.18,0.9, r'$P_z = 0$',fontsize=12)
fig.text(0.05,0.76, r'1NN',fontsize=12)
fig.text(0.05,0.50, r'2NN',fontsize=12)
fig.text(0.05,0.22, r'3NN',fontsize=12)
fig.text(0.07,0.86, r'(a)',fontsize=12)
fig.text(0.32,0.86, r'(b)',fontsize=12)
fig.text(0.575,0.86, r'(c)',fontsize=12)
fig.text(0.07,0.59, r'(d)',fontsize=12)
fig.text(0.32,0.59, r'(e)',fontsize=12)
fig.text(0.575,0.59, r'(f)',fontsize=12)
fig.text(0.07,0.31, r'(g)',fontsize=12)
fig.text(0.32,0.31, r'(h)',fontsize=12)
fig.text(0.575,0.31, r'(i)',fontsize=12)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
bar = fig.colorbar(l1, cax=cbar_ax, ticks=[-0.005, 0.000, 0.005, 0.010, 0.015])
bar.ax.set_yticklabels(['-0.005', '0.00', '0.005', '0.010', '0.015'], fontsize = 12) # vertically oriented colorbar
fig.text(0.99,0.54, r'OMA ($\mu_{B}$)',fontsize=12)
fig.suptitle('Effect of Perplexity (per = 20)', fontsize=16)
#plt.savefig("soap_10_unfilled.png", dpi=300, bbox_inches='tight')
plt.show()
def get_pca_X(X):
pca = decomposition.PCA(n_components=2)
X_pca = pca.fit_transform(X)
return X_pca
# First Nearest Neighbours
X_pca_Pz_0_1NN = get_pca_X(X_0_1NN)
X_pca_Pz_0p25_1NN = get_pca_X(X_0p25_1NN)
X_pca_Pz_0p5_1NN = get_pca_X(X_0p5_1NN)
X_pca_Pz_0p75_and_1_1NN = get_pca_X(X_last_two_1NN)
# (Second Nearest Neighbours
X_pca_Pz_0_2NN = get_pca_X(X_0_2NN)
X_pca_Pz_0p25_2NN = get_pca_X(X_0p25_2NN)
X_pca_Pz_0p5_2NN = get_pca_X(X_0p5_2NN)
X_pca_Pz_0p75_and_1_2NN = get_pca_X(X_last_two_2NN)
# Third Nearest Neighbours)
X_pca_Pz_0_3NN = get_pca_X(X_0_3NN)
X_pca_Pz_0p25_3NN = get_pca_X(X_0p25_3NN)
X_pca_Pz_0p5_3NN = get_pca_X(X_0p5_3NN)
X_pca_Pz_0p75_and_1_3NN = get_pca_X(X_last_two_3NN)
fig, axs = plt.subplots(nrows=3, ncols=3, figsize=(12, 10),subplot_kw={'xticks': [], 'yticks': []})
# set the colormap
cmap = plt.cm.get_cmap('viridis')
# Uncomment to try different colormaps. Note: Make sure to use perceptually uniform sequential colormaps
#cmap = plt.cm.get_cmap('plasma_r')
#cmap = plt.cm.get_cmap('inferno')
#cmap = plt.cm.get_cmap('magma')
# set the marker
marker = '^'
# First nearest-neighbours
l1 = axs[0][2].scatter(X_pca_Pz_0p75_and_1_1NN[:,0:1], X_pca_Pz_0p75_and_1_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_1NN)#, cmap = cmap,edgecolors='none')
l2 = axs[0][1].scatter(X_pca_Pz_0p5_1NN[:,0:1], X_pca_Pz_0p5_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_1NN)#, cmap=cmap,edgecolors='none')
l3 = axs[0][0].scatter(X_pca_Pz_0_1NN[:,0:1], X_pca_Pz_0_1NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_1NN)#, cmap=cmap,edgecolors='none')
l1.set_facecolor('none')
l2.set_facecolor('none')
l3.set_facecolor('none')
# Second nearest-neighbours
l4 = axs[1][2].scatter(X_pca_Pz_0p75_and_1_2NN[:,0:1], X_pca_Pz_0p75_and_1_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_2NN)#, cmap=cmap,edgecolors='none')
l5 = axs[1][1].scatter(X_pca_Pz_0p5_2NN[:,0:1], X_pca_Pz_0p5_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_2NN)#, cmap=cmap,edgecolors='none')
l6 = axs[1][0].scatter(X_pca_Pz_0_2NN[:,0:1], X_pca_Pz_0_2NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_2NN)#, cmap=cmap,edgecolors='none')
l4.set_facecolor('none')
l5.set_facecolor('none')
l6.set_facecolor('none')
# Third nearest-neighbours
l7 = axs[2][2].scatter(X_pca_Pz_0p75_and_1_3NN[:,0:1], X_pca_Pz_0p75_and_1_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_last_two_3NN)#, cmap=cmap,edgecolors='none')
l8 = axs[2][1].scatter(X_pca_Pz_0p5_3NN[:,0:1], X_pca_Pz_0p5_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0p5_3NN)#, cmap=cmap,edgecolors='none')
l9 = axs[2][0].scatter(X_pca_Pz_0_3NN[:,0:1], X_pca_Pz_0_3NN[:,1:2], s=100, alpha=0.9, marker = marker, c = Y_0_3NN)#, cmap=cmap,edgecolors='none')
l7.set_facecolor('none')
l8.set_facecolor('none')
l9.set_facecolor('none')
plt.subplots_adjust(left=0.1,
bottom=0.1,
right=0.9,
top=0.88,
wspace=0.25,
hspace=0.2)
fig.text(0.65,0.9, r'$P_z \geq 0.75$',fontsize=12)
fig.text(0.42,0.9, r'$P_z = 0.5 $',fontsize=12)
fig.text(0.18,0.9, r'$P_z = 0$',fontsize=12)
fig.text(0.05,0.76, r'1NN',fontsize=12)
fig.text(0.05,0.50, r'2NN',fontsize=12)
fig.text(0.05,0.22, r'3NN',fontsize=12)
fig.text(0.07,0.86, r'(a)',fontsize=12)
fig.text(0.32,0.86, r'(b)',fontsize=12)
fig.text(0.575,0.86, r'(c)',fontsize=12)
fig.text(0.07,0.59, r'(d)',fontsize=12)
fig.text(0.32,0.59, r'(e)',fontsize=12)
fig.text(0.575,0.59, r'(f)',fontsize=12)
fig.text(0.07,0.31, r'(g)',fontsize=12)
fig.text(0.32,0.31, r'(h)',fontsize=12)
fig.text(0.575,0.31, r'(i)',fontsize=12)
fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
bar = fig.colorbar(l1, cax=cbar_ax, ticks=[-0.005, 0.000, 0.005, 0.010, 0.015])
bar.ax.set_yticklabels(['-0.005', '0.00', '0.005', '0.010', '0.015'], fontsize = 12) # vertically oriented colorbar
fig.text(0.99,0.54, r'OMA ($\mu_{B}$)',fontsize=12)
#plt.savefig("soap_pca.png", dpi=300, bbox_inches='tight')
plt.show()